load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; radius to degree
r2d = 45.0/atan(1.0)

if (.not. isvar("configFile")) then
    configFile = "/public/users/wrf/Codes/WRF_scripts/WRF_Post/StationLists/wind_stnlist.csv"
end if

lines = readAsciiTable(configFile,1,"string",1)
delim = ","

stnnum    = ndtooned(str_get_field(lines,1,delim))
stnlist    = ndtooned(str_get_field(lines,2,delim))
stnlon    = ndtooned(tofloat(str_get_field(lines,3,delim)))
stnlat    = ndtooned(tofloat(str_get_field(lines,4,delim)))
stndomain = ndtooned(str_get_field(lines,5,delim))
nstns     = dimsizes(stnlist)

levels             = (/10, 30, 50, 70/)
;levels             = (/10, 30, 50, 70, 80, 120/)
nlevels    = dimsizes(levels)

;;;;;;;;;;;;;;;;;;;;;;
; Env and Settings   ;
;;;;;;;;;;;;;;;;;;;;;;


nlevel = dimsizes(levels)

if (.not. isvar("output_prefix ")) then
    output_prefix = "~/wrfpost/test/station_sswdn"
    if (.not. isfilepresent(output_prefix)) then
        system("mkdir " + output_prefix)
    end if
end if

if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012103120/wrfout_d02_2012-10-31_12:00:00"
end if

; default output time interval is 30min
if (.not. isvar("output_interval")) then
    print("output_interval is undefined")
    output_interval = 0.5 ; output 1h precip
end if

a = addfile(input_file + ".nc", "r")
filename = systemfunc("basename " + input_file)

if (.not. isvar("start_date")) then
    start_date = substring(filename, 11, 23)
end if

domain = substring(filename, 7, 9)
print(domain)

times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
ntimes = dimsizes(hours)
print(ntimes)

; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")

wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
dt = round(output_interval / wrf_dt, 3)

;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop for each station ;
;;;;;;;;;;;;;;;;;;;;;;;;;

do istn = 0, nstns - 1
    print(stndomain(istn))
    if (.not. isStrSubset(filename, stndomain(istn))) then
        continue
    end if


    output_filename = output_prefix + "/wind_output_" + stnlist(istn) +  "_" + domain +  "_" + start_date + ".txt"
    ;;; latlon to XY ;;;
    res = True
    res@returnInt = False   ; return real values
    loc = wrf_user_ll_to_ij(a, stnlon(istn), stnlat(istn), res) - 1

    sxc0 = floattointeger(loc(0))
    syc0 = floattointeger(loc(1))

    syc1 = syc0 + 1
    sxc1 = sxc0 + 1

    print ("X/Y location is: " + loc)
    print("X0 is: " + sxc0)
    print("Y0 is: " + syc0)
    print("X1 is: " + sxc1)
    print("Y1 is: " + syc1)

    line_array = new(ntimes(0)*nlevel, "string")
    i_line = 0
    ;    line_array(0) = "time,                   u10,       v10,       ulvl1,     vlvl1,    hlvl1,      ulvl2,     vlvl2,    hlvl2,      alfa_u,    alfa_u2,   alfa_v,    alfa_v2,   u70,       u70_2,     v70,       v70_2"
    ;    i_line = i_line + 1

    ;;; loop for time steps ;;;
    ;    do istep = 0, ntimes(0) - 1


    ;;; Load Terrain Data ;;;
    yy = (/syc0, syc1/)
    xx = (/sxc0, sxc1/)
    kk = (/0, 1, 2/)

    ;;; load variables
    ; read variables from wrfout ncfile
    ter         = wrf_dyy_getvar(a,"ter",0,kk,yy,xx)
    z           = wrf_dyy_getvar(a,"z",-1,kk,yy,xx)
    pres_hpa    = wrf_dyy_getvar(a, "pressure", -1, kk, yy, xx)
    uvmet = wrf_dyy_getvar(a, "uvmet", -1, kk, yy, xx)
    uvmet10 = wrf_dyy_getvar(a, "uvmet10", -1, kk, yy, xx)

    ; interpolate to the site location
    tmp_uv = rm_single_dims(linint2( xx, yy, uvmet, False, loc(0), loc(1), 0))
    tmp_uv10 = rm_single_dims(linint2( xx, yy, uvmet10, False, loc(0), loc(1), 0))
    tmp_ter = rm_single_dims(linint2( xx, yy, ter, False, loc(0), loc(1), 0))
    tmp_z = rm_single_dims(linint2( xx, yy, z, False, loc(0), loc(1), 0))
    tmp_pres = rm_single_dims(linint2( xx, yy, pres_hpa, False, loc(0), loc(1), 0))

    ; do the vertical interpolation
    uv_arr = wrf_uv_intrp_vert(tmp_ter, tmp_z, tmp_uv, tmp_uv10, levels)

    ; calculate wind speed and direction
    spd_arr = sqrt(uv_arr(0,:,:)^2 + uv_arr(1,:,:)^2)
    dir_arr = round(atan2(uv_arr(0,:,:), uv_arr(1,:,:)) * r2d + 180, 3)

    ;;; END OF WIND EXTRACTION ;;;

    ;;; write to txt file ;;;
    header = "ID,       TIME"

    ; station number and time
    printArray = stnnum(istn) + ", " + timeStr_cst
   
    ; wind speed and direction
    do ilvl = 0, nlevel - 1
        header = header + ", " + sprinti("%3i", levels(ilvl)) + "sp, " + sprinti("%3i", levels(ilvl)) + "dir"
        printArray = printArray + ", " + sprintf("%5.2f",spd_arr(:,ilvl)) + ", " + sprinti("%3i",dir_arr(:,ilvl))
    end do
    allPrintArray = new(ntimes + 1, "string")
    allPrintArray(0) = header
    allPrintArray(1:ntimes) = printArray

    asciiwrite(output_filename, allPrintArray)
    print("Creating: " + output_filename)

    delete(ter)
    delete(z)
    delete(pres_hpa)
    delete(printArray)
    delete(allPrintArray)
end do
